# Libraries
library(tidyverse)
library(sf)
library(USAboundaries)
library(ggthemes)
library(cowplot)
library(sp)
library(leaflet)
A minimum threshold of 10 distinct measurements with at least 1 measurement every 5 years was applied to filter for only the regularly monitored wells. investigate regularly monitored wells
# Data from USGS National Water Information System (NWIS)
# Arizona state shape (CRS = 5070)
az = us_states() %>%
filter(name == 'Arizona') %>%
st_transform(5070)
# add 'well #' for each unique well
usgs = az_nwis_unique_sites
usgs$wellid = paste("well", 1:nrow(usgs))
# remove wells located at same lat/long but have different well ID
usgs_spatial = usgs %>%
st_as_sf(coords = c('long_nad83', 'lat_nad83'), crs = 4269) %>%
st_transform(5070) %>%
as_Spatial() %>%
remove.duplicates() %>%
st_as_sf() %>%
group_by(wellid)
# match 'well #' from spatial data frame to corresponding well in time series data frame
usgs_time = az_nwis_all %>%
group_by(site_id)
usgs_time = left_join(usgs_time, select(usgs, site_id, wellid), by = "site_id")
usgs_time = usgs_time %>% filter(wellid %in% usgs_spatial$wellid)
# FILTER BY THRESHOLD - at least 10 measurements and at least 1 measurement every 5 years
thresh = usgs_time %>%
group_by(wellid) %>%
filter(measurement_dist >= 10) %>%
mutate(measure_period = year - lag(year)) %>%
filter(any(measure_period > 5))
usgs_time = usgs_time %>% filter(!wellid %in% thresh$wellid, measurement_dist >= 10) %>%
mutate(measure_period = year - lag(year))
usgs_spatial = usgs_spatial %>%
filter(!wellid %in% thresh$wellid) %>%
filter(measurement_dist >= 10)
# Data from Arizona Department of Water Resources (ADWR)
# add 'well #' for each unique well
state = adwr_unique_sites %>%
filter(date > 1960-01-01, date < as.Date.character('2000-01-01'))
state$wellid = paste("well", 1:nrow(state))
# remove wells located at same lat/long but have different well ID
state_spatial = state %>%
st_as_sf(coords = c('long_nad83', 'lat_nad83'), crs = 4269) %>%
st_transform(5070) %>%
as_Spatial() %>%
remove.duplicates() %>%
st_as_sf() %>%
group_by(site_id)
# match 'well #' from spatial data frame to corresponding well in time series data frame
state_time = adwr_all %>%
group_by(site_id)
state_time = left_join(state_time, select(state, site_id, wellid), by = "site_id")
state_time = state_time %>% filter(wellid %in% state_spatial$wellid)
# FILTER BY THRESHOLD - at least 10 measurements and at least 1 measurement every 5 years
thresh_state = state_time %>%
group_by(wellid) %>%
filter(measurement_dist >= 10) %>%
mutate(measure_period = year - lag(year)) %>%
filter(any(measure_period > 5))
state_time = state_time %>% filter(!wellid %in% thresh_state$wellid, measurement_dist >= 10) %>%
mutate(measure_period = year - lag(year))
state_spatial = state_spatial %>%
filter(!wellid %in% thresh_state$wellid) %>%
filter(measurement_dist >= 10)
I decided to include these increments of 100ft as it is visually easier to track individual wells



Located just north of Harquahala AMA boundary





Well 3450 was monitored continuously since ~1965 

Wells 1226 & 2710 had measurements recorded only for a couple years in the late 1990s


Northwest of Flagstaff and northeast of Kingman

shows a continual decline in DTW from ~1975 - 2000

Well 1411 is nearby to well 1669 and shows the same trend as 1669, 1411 has a max depth of 830 ft, just under 1000ft marker

Using the same threshold applied to the USGS data of 10 distinct measurements with at least 1 measurement every 5 years, the total number of State monitored wells was significantly reduced.
# PLOTS A BUFFER AROUND SPECIFIED WELL AND LOCATES OTHER WELLS INSIDE THE BUFFER AREA
buffer_fun = function(df, well, buff, state) {
buffer = st_buffer(df[well,], buff)
near1 = st_intersection(df[,], buffer) %>% filter(measurement_dist >= 10)
plot = ggplot() +
geom_sf(data = state) +
geom_sf(data = buffer, fill = NA) +
geom_sf(data = near1, col = "red", size = .5) +
labs(caption = paste(nrow(near1), 'wells')) +
theme_void() +
theme(plot.caption = element_text(size = 22, face = "bold", hjust = 0.5))
print(plot)
return(near1)
}
# RETURNS DATAFRAME OF WELLS IN DTW RANGE (MIN - MAX)
dtw_range = function(df, min, max) {
df = df %>% filter(dtw <= max, dtw >= min) %>%
mutate(sd = sd(dtw)) %>%
arrange(desc(date))
return(df)
}
# PLOTS AN INDIVIDUAL WELL BY WELL ID (DTW TIME SERIES)
plot_well = function(df_time, id) {
wells = df_time %>% filter(wellid == paste('well', id))
plot = ggplot(data = wells, aes(x = date, y = dtw_ft)) +
geom_line(aes(y = dtw_ft, col = wellid), size = 2) +
scale_y_reverse() +
labs(title = paste('Well', id),
caption = paste('Min depth:', min(wells$dtw_ft), '\nMax depth:',
max(wells$dtw_ft), '\nNumber of measurements:',
wells$measurement_dist),
x = 'Year',
y = 'DTW (ft)') +
theme_bw() +
theme(plot.title = element_text(face = 'bold',color = 'black', size = 18, hjust = 0.5),
axis.text.x = element_text(color="black", size=14),
axis.text.y = element_text(color="black", size=14),
axis.title.x = element_text(face = 'bold', color="black", size=16),
axis.title.y = element_text(face = 'bold', color="black", size=16),
plot.caption = element_text(face = 'bold', color = 'black', size = 14),
panel.grid.major = element_line(colour = "#808080"),
panel.grid.minor = element_line(colour = "#808080", size = 1),
legend.position = 'none')
print(plot)
}
# PLOTS MULTIPLE WELLS IN DTW RANGE
multi_well_plot = function(df, df_time, min, max) {
df = df %>% filter(dtw <= max, dtw >= min) %>%
arrange(desc(date))
df_time = df_time %>% filter(wellid %in% df$wellid)
plot1 = ggplot(data = df_time, aes(x = date, y = dtw_ft)) +
geom_line(aes(y = dtw_ft, col = wellid), size = 1) +
scale_y_reverse() +
labs(x = 'Year',
y = 'DTW (ft)') +
theme_bw() +
theme(plot.title = element_text(face = 'bold',color = 'black', size = 18, hjust = 0.5),
axis.text.x = element_text(color="black", size=14),
axis.text.y = element_text(color="black", size=14),
axis.title.x = element_text(face="bold", color="black", size=16),
axis.title.y = element_text(face="bold", color="black", size=16),
panel.grid.major = element_line(colour = "#808080"),
panel.grid.minor = element_line(colour = "#808080", size = 1),
legend.position = 'none')
print(plot1)
return(df)
}
# PLOTS WELLS DTW WITHIN A SPECIFIED RANGE, COLOR GRADIENT BY DEPTH, RETURNS CORRESPONDING DATA FRAME
map_dtw = function(state, df, min, max) {
df = df %>% filter(dtw <= max, dtw >= min) %>%
mutate(sd = sd(dtw)) %>%
arrange(dtw)
plot1 = ggplot() +
geom_sf(data = state, size = 1, col = 'black') +
geom_sf(data = df, aes(col = dtw), size = 1.5) +
theme_void() +
theme(plot.title = element_text(face = "bold", color = "black", size = 18),
plot.subtitle = element_text(size = 12),
plot.caption = element_text(face = 'bold', size = 12),
legend.title = element_text(color="black", size=14),
legend.text = element_text(color="black", size=14)) +
scale_colour_gradient() +
labs(title = paste('DTW', min, '-', max, 'ft'),
fill = 'Number of wells',
col = 'Depth to Water (ft)')
print(plot1)
return(df)
}